Counts

Poisson Distribution

The simplest way to model a count is with a Poisson distribution

Poisson Distribution

The expected value of \(Y|X\) is

\[ E(Y|X) = e ^{(b_0 + b_1X_1 +b_2X_2...)} \]

poisson_func <- function(x,beta0,beta1) {
  theta <- exp(beta0 + beta1*x)
}
curve(poisson_func(x,0,1),xlim=c(-1,5), lwd=2,ylab="poisson_func function", las=1)
curve(poisson_func(x,-2,1),xlim=c(-1,5), lwd=2, lty=2, col="red", add=TRUE)
curve(poisson_func(x,2,1),xlim=c(-1,5), lwd=2, lty=2, col="blue", add=TRUE)

legend(-1,150,c(expression(paste(beta[0]," = 0")),
                expression(paste(beta[0]," = -2")),
                expression(paste(beta[0]," = 2"))), 
       lwd=2, lty=1:3, col=c("black","red","blue"))

Poisson Distribution

poisson_func <- function(x,beta0,beta1) {
  theta <- exp(beta0 + beta1*x)
}
curve(poisson_func(x,3,-2),xlim=c(-1,5), lwd=2,ylab="poisson_func function", las=1)
curve(poisson_func(x,3,0),xlim=c(-1,5), lwd=2, lty=2, col="blue", add=TRUE)
curve(poisson_func(x,3,2),xlim=c(-1,5), lwd=2, lty=2, col="red", add=TRUE)

legend(3,150,c(expression(paste(beta[1]," = 0")),
               expression(paste(beta[1]," = -2")),
               expression(paste(beta[1]," = 2"))), 
       lwd=2, lty=1:3, col=c("blue","black", "red"))

Poisson Overdispersion

The poisson distribution assumes that the variance is equal to the mean:

set.seed(100)
x<-rpois(4000, 25)
var(x)

[1] 24.87314

mean(x)

[1] 25.01825

Poisson Overdispersion

The poisson distribution assumes that the variance is equal to the mean:

But this is rarely true with real world data:

  # count of terrorist attacks per country-year
mean(gtd_counts$nkill)

[1] 2.238636

var(gtd_counts$nkill)

[1] 63.0513

data.frame(dist = dist_poisson(lambda =mean(gtd_counts$nkill)))|>
  ggplot(aes(xdist = dist))+
  stat_slab(alpha=.7) +
  stat_slab(data=gtd_counts, aes(x=nkill), 
               fill='lightblue', 
               inherit.aes = FALSE) +
  theme_bw() +
  labs(title = 'Terrorist attacks in 2018 compared to a poisson distribution with the same mean') 

Poisson Overdispersion

• Models with more variance than expected are called “overdispersed”, and this causes us to underestimate our standard errors.

  • This can also happen with other distributions, but its probably less common because most distributions allow the variance to be different from the mean.

• This is may be the result of different data generating processes, i.e.: one pattern for states with terrorist insurgencies

  • In some cases there may be a way to fix this with a different set of parameters, but there’s no guarantee.

Options for handling overdispersion

Testing Overdispersion

The performance package has a function to check for overdispersion using simulated residuals from a count model. The null hypothesis here is no overdispersion:

library(performance)
check_overdispersion(model)

# Overdispersion test

       dispersion ratio =   18.614
  Pearson's Chi-Squared = 3201.581
                p-value =  < 0.001

We can also see this is addressed by the negative binomial model:

check_overdispersion(nb_model)

# Overdispersion test

 dispersion ratio = 0.689
          p-value = 0.864

Model fit

Since the only difference here is the addition of a single dispersion parameter, we can treat these two models as nested and use a likelihood ratio test to determine with the negative binomial model is a better fit. Unsurprisingly, it is:

pchisq(2 * (logLik(nb_model) - logLik(model)), df = 1, lower.tail = FALSE)

'log Lik.' 2.798821e-244 (df=5)

Poisson Prediction and interpretation

The negative binomial model and the poisson model both have the same inverse link function, so generating predictions here is the same:

• Calculate the linear prediction (plug in the values and multiple each by the coefficient)

• Exponentiate the result

  coefs<-coef(nb_model)
  exp(coefs['(Intercept)'] + 
        coefs['v2x_libdem'] * 0.74 + 
        coefs["any_attacksyes"] * 1 +
        coefs['e_gdppc'] * 59.6
  
      )

  (Intercept) 
     3.581757 

Count models are still nonlinear! So, like with logits, the marginal effect of X will depend in part on the other parameters.

Using a Random Effect

library(lme4)


mmodel<-glmer( nkill ~ 
                v2x_libdem + 
                any_attacks + 
                e_gdppc +
                (1|country_txt), data=gtd_counts, family='poisson')




performance::check_overdispersion(mmodel)

# Overdispersion test

       dispersion ratio = 0.012
  Pearson's Chi-Squared = 2.001
                p-value =     1